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Abstract: Biological systems are typically very complex and need to be reduced 
before they are amenable to a thorough analysis. Also, they often possess 
functionally important dynamic features like bistability. In model reduction, it 
is sometimes more desirable to preserve the dynamic features only than to recover 
a good quantitative approximation. We present an approach to reduce the order 
of a bistable dynamical system significantly while preserving bistability and the 
switching threshold. These properties are important for the operation of the system 
in the context of a larger network. As an application example, a bistable model 
for caspase activation in apoptosis is considered. 
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1. INTRODUCTION 

The importance of switch-like decisions in bi- 
ological processes has been revealed in a wide 
range of systems. Examples range from cell fate 
decisions via the MAP kinase cascade [Ferrell and 
Xiong, 2001] or apoptosis [Eifiing et at, 2004] to 
changes of metabolic parameters, as e.g. induced 
by the lac operon [Yildirim et ai, 2004]. Switch- 
like decisions are also a major feature of more 
complex dynamics, as found for oscillations in the 
cell cycle [Pomerening et ai, 2003]. Concerning 
the mathematical modeling, switch-like decisions 
are usually represented by bistable dynamical sys- 
tems. These systems have two stable steady states, 
and depending on initial conditions or external 
stimuli, converge to one of these steady states. 
Several theoretical approaches have been devel- 
oped to study the existence of two stable steady 
states in dynamical systems [Angeli and Sontag, 
2004, Eifiing et al, 2007]. 
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Another important issue in the modeling of bio- 
logical systems is model reduction [Conzelmann 
et al, 2004]. Due to the complexity of biologi- 
cal systems, typically only simplified models are 
amenable to a thorough computational analysis. 
Methods for model reduction which are theoret- 
ically founded provide important means to ap- 
proximate a detailed description of a system by 
a simpler model. 

Complex biological systems can often be viewed 
as a set of interconnected modules, each playing 
a specific role. In this case, two goals for model 
reduction of single modules can be distinguished. 
The first goal is to get a good quantitative approx- 
imation of the original model, such that a solution 
of the reduced model will differ as little as possible 
from a solution of the full model in the relevant 
variables (e.g. model outputs). For biochemical 
systems, this is often achieved via a time scale 
separation [e.g. Roussel and Eraser, 2001]. The 
second goal focusses on preserving the qualitative 
dynamics of a module and its role in the network 



[Dano et at, 2006]. In large biological networks, 
which typically are robust against fluctuations, 
the exact trajectories of the module might not 
be relevant, provided its qualitative behavior is 
preserved such that it can maintain its role in the 
network. When pursuing the latter goal, one can 
anticipate a much larger reduction than for the 
first one. For modules having the role of switches, 
the goal is thus to reduce the model while pre- 
serving bistability and the quantitative properties 
that are important for the module's operation, like 
the switching threshold. 

Based on a method introduced by Schmidt and Ja- 
cobsen [2004] , which allows to compute a measure 
of the contribution of each state variable in the 
model to bistability, we present an approach for 
model reduction preserving bistability in switch- 
like systems. The method is applied to a bistable 
system involved in programmed cell death. We 
show that the reduction preserves not only bista- 
bility. but also quantitative properties like the 
linear approximation to the manifold separating 
the two regions of attraction of the stable steady 
states, which is related to the switching threshold 
of the system. 



2. THE CASPASE ACTIVATION MODEL 

We study a model for caspase activation, a ma- 
jor part of apoptosis, developed by Eifiing et al. 
[2004]. Apoptosis, also called programmed cell 
death, is a signaling program that leads the cell to 
commit suicide under appropriate internal and/or 
external stimuli. It provides a living organism 
with means to remove infected, malfunctioning or 
simply unneeded cells to ensure its survival. Mal- 
function of apoptotic signaling has been detected 
in several diseases, including developmental de- 
fects, neurodegeneration and cancer [Danial and 
Korsmeyer, 2004]. Therefore it is also of medi- 
cal interest to understand the signaling network 
which regulates apoptosis. 

Cell death is a switch like decision: based on the 
input signal, the cell has to decide whether to 
undergo apoptosis or to stay alive. There is no 
gradual response. The role of the switching ele- 
ment in apoptosis is taken by the caspase cas- 
cade. One distinguishes initiator caspases which 
receive the stimuli and effector caspases actually 
carrying out apoptosis. Bistability in the model 
arises from a positive feedback loop between the 
effector caspases and the initiator caspases, in 
connection with the specific inhibitors for each 
type of caspases. 

The model as developed by Eifiing et al. [2004] 
already contains several simplifications, in partic- 
ular several types of initiator and effector caspases 




Fig. 1. Representation of the caspase activa- 
tion model. Constitutive degradation reac- 
tions are not displayed. 

are combined in one species each, and the same 
applies to several types of inhibitors of the effector 
caspases. Furthermore, the external stimulus is 
not explicitly included as an input, but instead 
the initial amount of activated initiator caspases 
resulting from the stimulation is considered as 
input to the model. 

The species present in the model are the initiator 
caspase 8 and the effector caspase 3, both in active 
(C8a, C3a) and inactive (C8, C3) forms, and the 
inhibitors lAP and CARP as well as their complexes 
with active caspase 3 and 8, respectively. Figure 1 
displays the species involved in the model and 
the reactions among them. The variables in the 
mathematical model denote the molecule numbers 
per cell of the various species (Table 1). 

Xi = [C8] xr, = [lAP] 

X2 = [C8a] x-g = [C3a_IAP] 

xz = [03] X7 = [CARP] 

X4 = [C3a] X8 = [C8a_CARP] 

Table 1. Model state variables. 

The equations for the caspase activation model are 
given in Table 2 with parameter values as shown 
in Table 3 [Eifiing et al., 2004]. These parameter 
values have been collected from literature. For 
simplicity, we consider state variables (given in 
molecules per cell) as dimensionless. 

The model yields two stable steady states, which 
can be identified with the living state, where 
the concentrations of all active caspases — both 
free and bound — are equal to zero, and with 

XI = —k2XlX4 — kgXl + fc_9 

X2 = k2X4Xi — k5X2 — feiia;2a;7 + k-iixs 
X3 = —kiX2X3 — kioxa + fc_io 

X4 = klX2X3 — k3X4X5 + k-3Xe — k6X4 

X5 = —k3X4XS + k-3X6 — k4X4X5 — ksXS + k-8 

xe = k3X4X5 - k-3XQ - k-jXQ 

±7 = -kllX2X7 + k-llXs - k\2X7 + fc-12 

is = kiiX2X7 — k-iixs — kisxs 

Table 2. Model equations 
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fcg 3.9 • 10-3 

fcio 3.9 ■ 10-3 

Table 3. Parameter values [min-^] 

the apoptotic state where an almost complete 
activation of the caspases determines the cell 
to undergo apoptosis. The model captures the 
switch-like behavior occurring during apoptosis 
very well and is thus suitable to our needs. 

Furthermore, we have a threshold manifold sep- 
arating the regions of attraction, as well as an 
unstable steady state on the threshold manifold 
which we refer to as the decision state, as it is 
relevant in the decision about the fate of the cell. 



3. DETERMINATION OF RELEVANCE 

When determining the relevance of each state 
variable in the system to bistability, the unstable 
steady state, or decision state, plays a crucial 
role. The boundaries of the bistability region 
in parameter space arc marked by bifurcations 
of the unstable steady state and thus can in 
principle be detected by a classical bifurcation 
analysis. Although bifurcation analysis reveals the 
influence of parameters on the bistability, it does 
not give the relevance of each variable and how 
these are interconnected to generate bistability. 

Based on a linear approximation of the model 
around the decision state, Schmidt and Jacobsen 
[2004] perturb the influence of one state variable 
on the other variables of the system. In the un- 
perturbed case, the linear approximation is un- 
stable, since the decision state is unstable. One 
now searches for a perturbation that renders the 
decision state stable, which would roughly cor- 
respond to reaching a bifurcation point in the 
original nonlinear system. The magnitude of a 
perturbation found in this way is then a measure 
of the relevance to bistability of the perturbed 
variable with its connections to other variables. 

To make this approach precise, consider a system 
given by the differential equation 

|. = /W. (2) 

with X G R". We assume that the system is 
bistable and has an unstable steady state Xd, the 
decision state for the bistable behavior. Setting 
Ax = x — Xd, the linearization around the decision 
state is given by 



Fig. 2. Closed feedback loop with one feedback 
path perturbed 

j^{Ax)=AAx, (3) 

with A = §|(a;d). Note that, by our assumption, 
the system (3) is unstable. 

The linearized system is then considered 
closed feedback loop, in the sense that all intercon- 
nections between the state variables are put into 
a virtual feedback path. Breaking this feedback 
path yields the control system 

L: ^(Ax) = AAx + (A- A)u, 
at 

with u G M" and 

fan \ 

A = 

i.e. A contains only the elements on the diagonal 
of A. By setting u = Ax, the feedback path is 
closed again and one obtains the original linear 
system (3). 

Following Schmidt and Jacobsen, we assume that 
the matrix A is stable, i.e. an < iov i = 1, . . . ,n. 

This implies that the instability of the closed 
loop system (3) is due to interconnections among 
variables. 

The computation is then done as follows: for each 

state variable Xi, a perturbation is introduced 
into the feedback path of this variable (Fig. 2). 
This yields the system 

Li{ei) : ^(Ax) = AAx + {Ai - Ai)€iAxi, 

where Ai and A^ are the i-th column of A and A, 
respectively. 

We are now searching for the minimal perturba- 
tion that will stabilize the system Lj and define 
the value as 

Ci = minjcj > | Li(ej) or Li{—ei) is stable}. 

(4) 

If the minimum does not exist, set 6^ = 00. The 
higher the value of ej, the more difficult it becomes 
to perturb the connections among the considered 
state variable Xi and the remaining variables such 
that instability is lost, and the less relevant the 



• CSa > > > > > C8SICSRF 

in' - • f'3a C3a_IAP • 

• • 



> 

10° 1AP ^ABP , 

C3 • • ; 

1°"' - C8 - 
I • ^ ^ ^ ^ ^ ^ ^ : 

1 2 3 4 5 6 7 8 

Components of caspase activation modei 

Fig. 3. Relevance of components to bistability in 
the caspase activation model 

variable Xi is to bistability. Formally, we use the 
following definition of relevance. 

Definition 1. The relevance Ri of the state vari- 
able Xi to the bistability of the system (2) is 

-Ri = — , 

where ti is given by equation (4) . 

This computation has been implemented numer- 
ically in the Systems Biology toolbox for Matlab 
[Schmidt and Jirstrand, 2006] and has been ap- 
plied to the model of caspase activation described 
in Section 2. The result is shown in Fig. 3, and 
one distinguishes easily the relevant components 
C8a, C3a, C3a_IAP and C8a_CARP from the oth- 
ers, which are much less relevant to bistability. 

It might be possible to develop a formal ap- 
proach to compute e in equation (4) based on 
the Kharitonov theorem and its extensions [Bhat- 
tacharyya et ai, 1995]. However, the cost of com- 
puting a stabilizing perturbation would be similar 
to the one of a direct numerical search, since we 
consider only one perturbation at a time. Thus 
this idea is not pursued here. 

4. BISTABILITY PRESERVING MODEL 
REDUCTION 

Based on the measure of relevance presented in 
the previous section, we now develop a method 
of model reduction which preserves bistability in 
the reduced model. This reduction method is then 
applied to the caspase activation model (1). 

4-1 Method of model reduction 

The basic idea of the model reduction already 
used by Schmidt and Jacobsen [2004] is to retain 



only the state variables which have been identified 
as relevant to the bistability, and to neglect the 
dynamics of the other variables. In contrast to 
Schmidt and Jacobsen, who replace the neglected 
state variables in the remaining equations with 
their steady state values, we use their steady state 
map which is computed from the full model. Our 
approach will be justified in Section 4.2. 

The relevant state variables can be identified by 
choosing appropriate thresholds r > for the 
relevance and a > 1 which determines the gap 
size. 

Definition 2. The state variable Xi, i — 1, . . . ,n, 
is said to be relevant (resp., irrelevant) to bista- 
bility, if there exist r > and a > 1 such that 
Ri > ar (resp., Ri < ^r). 

Using this definition, the state x of the system (2) 
is subdivided as 

X ^ ixB.,Xl)'^, 

where Xji contains the relevant variables and xj 
the irrelevant ones. Note that a should not be 
chosen to small to get a clear distinction between 
relevant and irrelevant state variables. 

We then proceed as follows: 

(1) For each variable which is not relevant to 
bistability, compute its steady state map 
from the full model, i.e. solve the equation 

0^f{xH,xi) (5) 

for xj to find the steady state map 

xi=g{xR). (6) 

Computationally, this is similar to the di- 
vision of a system in fast and slow subsys- 
tems using a quasi-steady-state approxima- 
tion [Schauer and Heinrich, 1983], though we 
do not have fast and slow subsystems here, 
but rather subsystems that are relevant or 
irrelevant to bistability. The map / has to 
be invertible with respect to xj, which is the 
case in our example. Local invertibility can 
be checked via the implicit function theorem, 
but easily checkable conditions for global in- 
vertibility are not available in general. 

(2) Drop the differential equations for xi from 
the system and replace the components of xi 
in fn by their steady state map g(xR). 

The reduced model is thus given by 

xr = fR{xR,g{xR)). 

To apply the described method to the model of 
caspase activation (1), we choose r = 10 and 
a = 10. This yields the separation 

XR = {x2,X4,xg,xs)^ and xi = {xi,X3,X5,X7)'^ . 



For the first step in the reduction, one gets the 
steady state map for xi as 

k-9 k-io 



Xi = 



X5 
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(fcs + k4)X4 + ks 



X3 
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kl2 + kiiX2 



(7) 

and the dynamics for the reduced model are thus 
given by 

k-gk2X4 , {k-12 + k-llX8)kllX2 , , 
X2 = — ; ksX2 ; ; h K-ll^s 



X4, 
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(fc-12 + fc-iij:8)fciia;2 
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■ k-UXg — kl3X8- 



(8) 



Note that, in addition to the reduction in system 
dimension, also the dimension of the parameter 
space can be reduced: setting kg — kg/k2 and 
^10 — kio/ki, the number of independent param- 
eters in the model is reduced by two. 



4-2 Results of reduction and interpretation 

Analyzing the dynamics of the reduced model, one 
gets the following important results: 

• All steady states of the full model are repro- 
duced as a projection in the reduced model. 

• The steady states of the reduced model have 
the same local stability behavior, in particu- 
lar the same number of eigenvalues with pos- 
itive real parts as their corresponding steady 
states in the full model after linearization. 

• The linear approximation to the stable man- 
ifold at the decision state is the same up to 
projection for both the full and the reduced 
model. 

The first property is obtained automatically by 
keeping the steady state map (7) of the ne- 
glected variables, under the assumption that no 
two steady states are projected to the same state 
by the reduction. Then one directly achieves the 
preservation of all steady states, which is a nec- 
essary condition for preservation of bistability. At 
this point, the change we made to the original 
approach of Schmidt and Jacobsen [2004] is im- 
portant. Replacing the neglected variables with 
their constant steady state values does in general 
not preserve steady states other than the unstable 
decision state, in particular it does not preserve 
the two stable steady states in the caspase ac- 
tivation model. But to preserve bistability, it is 
not only necessary to preserve instability of the 
decision state, but one also needs to preserve the 
two stable steady states. The use of the steady 



state map provides a general way to preserve all 
steady states in the reduced model, which can- 
not be achieved by using constant values for the 
irrelevant variables. 

The second and third property depend on a clear 

separation in relevant and irrelevant state vari- 
ables and thus on the choice of the parameters 
r and a. Preservation of bistability is actually 
due to the second property: it guarantees that 
the living steady state and the apoptotic steady 
state are stable in the reduced model, while the 
third steady state, which is part of the threshold 
manifold, remains unstable. 

The third property gives some quantitative mea- 
sures that are reproduced exactly in the reduced 
model. The stable manifold of the decision state 
represents the switching threshold surface for the 
bistable system. Their linear approximations are 
equivalent in the original and the reduced model. 
This implies that we have not only recovered the 
qualitative trait of the system being bistable, but 
also quantitative measures like threshold values 
have been preserved when staying close to the 
decision state. It has to be expected though that 
the models will differ more when further away 
from the decision state, since the reduction was 
based on a local analysis at the decision state. 

Biologically, it is interesting that the method re- 
vealed the active forms of the caspascs as the 
relevant elements in generating bistability, and 
dropped all inactive forms and free inhibitors. 
Using the steady state maps for these allows the 
reduced system to retain its role as a biological 
switch in the apoptotic network. In contrast to 
a time scale separation, the removed elements do 
not show a much faster dynamics here. Yet their 
dynamics are not contributing to bistability, and 
thus can safely be discarded. The removed vari- 
ables are then considered as a pool of substances 
being in quasi steady state. 

A comparison of the full and the reduced model by 

numerical simulation yields the results displayed 
in Fig. 4, obtained by starting with different 
initial conditions of free C8a. This figure illustrates 
that the linear approximation to the threshold 
manifold (LATM) is the same up to projection 
for both models, as well as the unstable direction 
of the decision state. Furthermore, the trajectories 
are very similar when starting close to the decision 
state. 

However, we also observe that the reduced system 
has different dynamics when considered further 
away from the decision state. In particular the 
threshold in initial C8a required to undergo apop- 
tosis is lower in the reduced model. Note that 
for the full model, the linear threshold displayed 
in Fig. 4 is a very good approximation of the 
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Fig. 4. Trajectories of the full and the reduced 
model projected to a plane spanned by vec- 
tors pi and p2- These were chosen such that 
(a) the LATM gets projected to a line and (b) 
they represent approximately the amount of 
total C3a and C8a. respectively. 

threshold manifold, while the approximation is 
not that good for the reduced model. This is likely 

due to the higher nonlincaritics in the reduced 
model introduced via the steady state maps of the 
neglected components. 

Although the trajectories of the full and the 

reduced model are different when not starting 
close to one of the steady states, the results are 
still meeting our goals, since the focus of our study 
is more on preserving bistability than recovering a 
quantitatively similar model. In this way, we have 
identified a minimal set of components responsible 
for switch-like behavior. 

5. CONCLUSIONS 

We have studied a mathematical model of a cas- 
pase activation system involved in the initiation 
of apoptosis. The model is bistable to incorpo- 
rate the switch like behavior observed in the real 
system. We have computed the relevance of each 
state variable to bistability and observed that 
only four out of eight variables have a high rele- 
vance. Using this result, the model can be reduced 



significantly by retaining only the relevant state 
variables in the model equations plus the steady 
state map of the irrelevant variables. The reduced 
model retains the bistability of the original model 
and also quantitative features such as trajecto- 
ries close to the steady states and the linear ap- 
proximation to the threshold manifold separating 
the regions of attraction of the two stable steady 
states. 
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